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We present an effective hydrodynamic theory of electronic transport in graphene in the interaction- 
dominated regime. We derive the emergent hydrodynamic description from the microscopic Boltz¬ 
mann kinetic equation taking into account dissipation due to Coulomb interaction and find the 
viscosity of Dirac fermions in graphene for arbitrary densities. The viscous terms have a dramatic 
effect on transport coefhcients in clean samples at high temperatures. Within linear response, we 
show that viscosity manifests itself in the nonlocal conductivity as well as dispersion of hydrody¬ 
namic plasmons. Beyond linear response, we apply the derived nonlinear hydrodynamics to the 
problem of hot spot relaxation in graphene. 


Physics at long time and length scales can be conve¬ 
niently described within the hydrodynamic approach^. 
The appeal of this approach is hinged on its ability to de¬ 
scribe a wide range of physical system^^ using the same, 
relatively small set of quantities and equations governing 
their behavior. At the same time, the final form of the 
hydrodynamic equations varies from system to systenP^ 
reflecting the particular symmetries and other physical 
features of the problem. 

Traditional hydrodynamicJi^ describes the system in 
terms of the velocity field v. The equations describing 
the velocity field (e.g., the Euler equation in the case of 
the ideal liquid or the Navier-Stokes equation if dissipa¬ 
tion is taken into account) can be either inferred from 
symmetry arguments or derived from the Boltzmann ki¬ 
netic equation. Both approached require one to express 
the fluxes of conserved quantities (energy, momentum, 
etc.) in terms of v. In particular, the viscous terms ap¬ 
pearing in the Navier-Stokes equation can be traced to a 
particular approximation for the momentum flux (or the 
stress tensor) Wap- The specific form of Wap depends 
on whether one discusses a usual, Galilean-invariant or a 
relativistic, Lorentz-invariant system. 

Low-energy excitations in graphen^ present a most 
interesting case of a system that is neither Galilean- nor 
Lorentz-invariant. This poses a significant challenge in 
establishing the hydrodynamic description in grap hene, 
which has to be derived from first principlesSHIS. The 
resulting equations should account for physical processes 
at time and length scales that are much longer than the 
scales related to the microscopic processes responsible for 
equilibration of the system. The issue of scale separation 
is especially important in the vicinity of charge neutral¬ 
ity in clean graphene. Without interaction there is no 
“intrinsic” energy scale other than temperature. 

The interest in hydrodynamics in graphene has been 
underpinned by the tremendous promise for potential ap¬ 


plications, e.g., for optoelectronic^^oiau^ where the hydro- 
dynamic approach is particularly suit able for describing 
the low-frequency optical responsd^l22l Linearized hydro- 
dynamic equations provide effective tools for evaluating 
transport coefficien ts in gr aphene and graphene-based 
double-layer devicesl^^^. At the same time, novel ex¬ 
perimental technique^^IIIIHlll bring the studies of nonlin¬ 
ear effects an nonlocal transport phenomena in graphene 
is within reach, while improved fabrication methods have 
yielded ultra-clean sample^^. For example, graphene on 
hexagonal boron nitride has been shown to support as¬ 
tonishingly homogeneous charge densitie^^. 

In this paper we derive a hydrodynamic description 
of electronic transport in graphene in the collision- 
dominated regime, where the shortest time scale in the 
problem is provided by electron-electron interaction. On 
the contrary, time scales associated with potential disor¬ 
der are assumed to be the longest in the system. Conse¬ 
quently, disorder plays no role in our theory. Our deriva¬ 
tion is based on the quantum kinetic equation (QKE) 
approach, which has been previou^ used to derive the 
macroscopic linear response theoryl^. 

The transition from the microscopic, kinetic descrip¬ 
tion to the macroscopic, hydrodynamic equations is sim¬ 
plified by the so-called “collinear sc attering singularity” 
of the collision integraP B^ I ^s i ^^ i a lHgi] QKE, i.e. 

the observation that kinematic properties of the Dirac 
quasiparticles lead to a divergence in the collision in¬ 
tegral for scattering processes involving quasiparticles 
moving along the same d irection . Dynamical screening 
regularizes the divergenc d^^ * ^^ * ^^ -*, such that the result¬ 
ing generic relaxation rates in graphene contain a large 
factor T~^ oc |lnag| ^ 1, where ag = e^jeVg is the ef¬ 
fective coupling constant (here e is the effective dielec¬ 
tric constant of the substrate and Vg is the “speed of 
light” in graphene). Dependin g on th e substrate, the 
coupling constant may be smalP^EMol^ a„ < 1. There 
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are, however, three macroscopic current^Slls] that are not 
relaxed at times of order Tg-. (i) the energy current j 
(ii) the electric current j ; and (iii) the so-called imbalance 
currentP^ jj. 

The energy current j^ in graphene is equivalent to the 
total momentum of electrons and thus cannot be relaxed 
by electron-electron interaction. The electric current in 
graphene is determined by the velocity rather than the 
momentum and therefore is not a conserved quantity. 
However, it is conserved in the collinear scattering pro¬ 
cesses and hence the corresponding relaxation rate does 
not contain the logarithmic enhancement. Finally, the 
imbalance current j is proportional to the sign of the 
quasiparticle energy and to the velocity. Similarly to 
the electric current, it does not experience logarithmi¬ 
cally enhanced relaxation. The imbalance current is re¬ 
lated to the quasiparticle number or imbalance densitj^^^, 
n/ = n+ + ri-, where n+ and n_ are the particle num¬ 
bers in the upper (conduction) and lower (valence) bands. 
Neglecting the Auger processes, quasiparticle recombina¬ 
tion due to e.g. electron-phonon interaction, and three- 
particle collisions due to weak coupling, one hnds that 
n+ and n_ are conserved independently. In this case, 
which will be considered in the rest of the paper, not 
only the total charge density n = n+ — n_, but also the 
quasiparticle density nj is conserved. 

At times longer than Tg, physical observables can be 
described within the macroscopic - or hydrodynamic - 
approach. The existence of the three slow-relaxing modes 
in graphene implies a peculiar two-step thermalization. 

Short-time electron-electron scattering (at time scales 
up to Tg) establishes the so-called “unidirectional therma- 
lization’^l^ the collinear scattering singularity implies 
that the electron-electron interaction is more effective 
along the same direction. Within linear response,li^ one 
can express the non-equilibrium distribution function in 
terms of the three macroscopic currents j, j^, and jj. 
The currents can then be found from the macroscopic 
equations. The currents j and jj are not conserved and 
can be relaxed by the electron-electron interaction. Close 
to charge neutralit y, th e corresponding relaxation rates 
can be estimated ^ a^T ^ r“^. These rates 

enter the macroscopic equations as friction-like terms. 
The macroscopic linear response theory has the same 
form on time scales shorter or longer than Tee- 

Beyond linear response, the scattering processes char¬ 
acterized by the time scale Tee play an important role 
in thermalizing quasiparticles moving in different direc¬ 
tions and thus lead to establishing the local equilibrium. 
This is the starting point for derivation of the nonlin¬ 
ear hydrodynamics, which is valid at time scales much 
longer than Tee- In view of conservation of the particle 
number, energy, and momentum, as well as independent 
conservation of the number of particles in the two bands 
in graphene, we may w rite the local equilibrium distri¬ 
bution function aJ^^Ell 

fx.lir) = {1 + exp [/3(r)(eA,fc - ma(t-) - u{r)-k)]}~\ 

( 1 ) 


where = Xvgk denotes the energies of the electronic 
states with the momentum k in the band A = ±1, 
the local chemical potential, the local temperature is en¬ 
coded in /3(r) = 1/T{r), and u{r) is the hydrodynamic 
velocity field which we define below (this field should 
not be confused with quasiparticle velocities v). The 
distribution function 0 follows from the standard argu¬ 
ment similar to the Boltzmann’s H-theorercPl the equilib¬ 
rium state is characterized by time-independent entropy. 
The particular form Q takes into account the symmetry 
properties of the two-body electron-electron interaction 
and is valid for arbitrary single-particle spectrum. The 
latter means that Eq. Q relies on neither Galilean nor 
Lorentz invariance. 

Expanding the local equilibrium distribution function 
Q up to the leading order in deviations from the uni¬ 
form, equilibrium Eermi distribution, we recover the dis¬ 
tribution function used in the linear response theorj®. 
As we have already mentioned, this linearized distribu¬ 
tion has the same form also on time scales shorter than 
Tee- This is a property of the linear approximation. 
Should we attempt to find the subleading nonlinear terms 
in the distribution function for t < Tee, the result would 
not correspond to the Taylor expansion of Eq. Q. 

Assuming the local equilibrium Q for times t 3> Tee, 
we derive the nonlinear hydrodynamics in graphe ne sim i- 
larly to the standard Chapman-Enskog procedur^^EIHlll 
The important feature of our theory is the larger than 
usual number of hydrodynamic modes (densities of con¬ 
served quantities): total charge, energy, and quasipar¬ 
ticle imbalance densities and the energy current. The 
independence of these modes can be traced to the spe¬ 
cific feature of the quasiparticle spectrum in graphene: 
the inequivalence of velocity and momentum. 

Having derived the hydrodynamic equations, we turn 
to consider a representative example of nonlinear physics 
in graphene, the relaxation of a hot spot. By this we 
mean a particular non-equilibrium state of the system 
that is characterized by a locally elevated energy density. 
Such a state can be prepared with the help of a local 
probe or focused laser radiation. Evolving the system 
according to the hydrodynamic theor y, we find a rather 
surprising result. Although as expectecP^l^, the hot spot 
emits plasmonic waves that carry energy away, a nonzero 
excess energy density remains at the hot spot. Physically, 
this effect appears due to compensation between the pres¬ 
sure and the self-consistent electric (Vlasov) field, which 
leads to a quasi-equilibrium. Taking into account the 
dissipation leads to the decay of the quasi-equilibrium 
energy density at the hot spot. This decay however, is 
characterized by a longer time scale compared to the ini¬ 
tial emission of the plasmonic waves. At the same time, 
viscous effects lead to damping of the plasmonic waves 
themselves. 

The remainder of the paper is organized as follows. 
In Sec. m we develop the nonlinear hydrodynamic the¬ 
ory including dissipative terms starting from the QKE. 
In Sec. [H] we briefly discuss linear response in graphene. 
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Finally, Sec. |III| is devoted to nonlinear hydrodynamics 
in graphene. Here we present results on the relaxation 
dynamics of a hot spot obtained by a numerical integra¬ 
tion of the hydrodynamic equations. Technical details, 
e.g., the calculation of scattering rates for the dissipative 
terms, are relegated to appendices. 

I. HYDRODYNAMIC THEORY IN GRAPHENE 

In this Section we develop a hydrodynamic theory of 
transport in graphene in the collision-dominated regime. 
We begin with a short overview of the microscopic mech¬ 
anisms responsible for establishing the hydrodynamic 
regime. The resulting hydrodynamic equations are sum¬ 
marized in Sec. II PI 

A. From the microscopic theory to hydrodynamics 


(where in the absence of magnetic held the steady state 
can exist without disorder), electron-electron interaction 
alone is insufficient for this task. Similarly, in this paper 
we keep in mind that infrared divergencies should be cut 
by disorder. However, for physical observables, e.g., op¬ 
tical response, in the frequency window ^ w 3> 
the impurity scattering is irrelevant. 

We assume that local equilibrium is established at 
time scales of the order of Tee, i.e. the longest time 
scale associated with two-particle electron-electron inter¬ 
action. The corresponding length scale, Zhydro VgTge, 
dehnes the size of the local huid elemeniP. Note, that 

^hydro ^ l/{alT) » 1/T. 

Following the standard line of argumeniP, small devi¬ 
ations from the local equilibrium can be accounted for 
by introducing a small correction 6f to the distribution 
function Q 

/ = + Sf. 


1. Microscopic description 

Microscopically, the electronic system is governed by 
the Boltzmann kinetic equation 

/:/ = 5tee[/]-T-dl.h/-(/)^), (2) 

with the standard Liouvillian form in the left-hand side, 

C = dt-\- V ■ Vr -I- [eE -\- e{v x B)] ■ V^, (3) 

and the collision integral in the right-hand side. Scat¬ 
tering off potential disorder is described within the usual 
T-approximation with r^is being the disorder mean free 
time. Electron-electron interaction is described by the 
collision integral Stee[f]- 

In graphene, the electronic states can be labeled by the 
momentum k and the band index A = ±. These states 
are characterized by the energies exk = Xvgk and veloc¬ 
ities V = Xvgk/k, where Vg ~ 10®m/s. Hereafter we will 
work with the units where Vg = 1. Consequently, the 
distribution function can be denoted as / = fxkir). The 
angular average in the disorder part of the collision inte¬ 
gral is defined as the average over the direction of fc, 

-I-IT 

{f)^ = I f = l^hk- ( 4 ) 

^ —TT 

In the interaction-dominated regime, the scattering 
time due to electron-electron interaction is much smaller 
than the disorder scattering time 

7'ee ^ "Fiis- 

The same condition was previously used in the derivation 
of the linear response theory in graphene. Within linear 
response, the role of disorder is to establish the steady 
state. With the exception of the charge neutrality point 


Kinematic restrictions imposed by the linear spectr um in 
graphene lead to the collinear scattering singularitj^^fiEil 
in Stee [f] ■ While the singularity is regularized by screen¬ 
ing, most eigenmodes of Stee[f] decay at the shortest 
time scales Tn ~ Tee/I Inaol. As a result, within the lead- 
ing logarithmic approximatiorP*^ only three modes con¬ 
tribute to the hydrodynamics and we can parametrize Sf 
as 


where 


Sf = T 



Sfii) + 


Sfd) 



J=i 


(5a) 


(5b) 


5/(2) = 


VaV/3 

J^2 




and the three modes are 


(Si = 1, 412 = A, 43 = £/T. 


(5c) 


(5d) 


The non-equilibrium corrections to the distri¬ 
bution function should leave the conserved quantities 
unchangecP. As a result, the coefficient /i(^) = 0 (which 
could be understood as a shift of the velocity u), while 
the tensors g^g have to be traceless [otherwise the three 
terms in Eq. (5c) would shift the particle number density 


n, imbalance density n/, and energy density he, respec¬ 
tively] . 

The coefficients ha^ and are determined from the 
QKlPlil, which becomes a matrix equation in the re¬ 
stricted subspace of modes 4j- In what follows, we will 
use a short-hand notation 


SUelf] « -CSf, 


( 6 ) 
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where the matrix C corresponds to the linearized colli¬ 
sion integral. The technicalities of inverting the matrix 
C are discussed in Appendix]^ where we also relate the 
matrix collision integral to the diagrammatic calculation 
of conductivity and viscosity based on the Kubo formula. 

Finally, macroscopic equations describing electronic 
transport in graphene are obtained by integrating the 
kinetic equation with the distribution function ([^. In 
the hydrodynamic regime, i.e. at time scales much longer 
than Tee, the natural macroscopic variables are the modes 
that are not relaxed by electron-electron interaction. All 
non-conserved quantities should be expressed in terms of 
such “hydrodynamic” modes. In graphene, these include 
the densities n, ni, and ue, and the energy current Je- 
The electric and imbalance currents can then be found us¬ 
ing the equations of state. The emerging hydrodynamics 
is valid as long as the macroscopic quantities vary slowly 
on the scale /hydro set by interactions. 


2. Macroscopic quantities 

Most two-body electron-electron collisions in graphene 
leave the particle number in each band unchanged. This 
is the consequence of the linear dispersion relation. The 
only exception is given by the so-called Auger processes, 
where the direction of the momentum of all initial and 
final states in each scattering event is the same (i.e., all 
four states belong to the same straight line on the dis¬ 
persion cone). In the absence of disorder, the probabil¬ 
ity of Auger processes vanishes within the random phase 
approximation. Even if impurity-assisted processes are 
taken into account, the recombination rate due to Auger 
processes remains small. Other processes that may con¬ 
tribute to quasiparticle recombination include electron- 
phonon interaction (by means of either two-phonon or 
impurity assisted scattering) and three-particle collisions. 
All these processes introduce parametrically small relax¬ 
ation rateJi^ (close to charge neutrality, at least of order 
QfgT). Here we will neglect recombination and assume 
the densities n± to be conserved independently. 

The particle and energy densities n± and ue can be 
calculated with the help of the distribution function in a 
standard way 


where A^ = 4 accounts for the spin and valley degeneracy 
in graphene. In Eq. (7c) we measure the energy density 


He with respect to heo, 


n-Eo — / £-.fc 

Jk 


,kj 


lk<A 


( 8 ) 


which is the total energy density at charge neutrality 
and zero temperature. The ultra-violet cut-off A must 
be formally included [also in Eq. ([7c|)]. However, it drops 
out of the physical results. 

The densities of the conduction and valence bands, 


Eq. (7a) and Eq. (7bI, can be combined into the total 


charge and imbalance densities 


n = n_|_ — n_, 

nj = n+ -I- n_. 

The macroscopic currents are defined 


3+ — / '^+.,kf+.kj 
Jk 


(9a) 

(9b) 

(lOa) 


3- = / 
Jk 


Je — 


/ £\,kVx,kf\,k- 

J\M 


(10b) 


(10c) 


The electron and hole currents, Eqs. (10a) and (10b) can 


be combined into the electric and imbalance (or quasi¬ 
particle) currents 


3 =3+-3-, (lOd) 

3i=j++3-, (lOe) 

In graphene the energy current j e is equivalent to the 
momentum and is conserved, while the electric and im¬ 
balance currents can be damped by electron-electron in¬ 
teraction. 


n+ 



(7a) 


B. Generalized Euler equation 


n_ 


/ (i-/-,fe), 

Jk 


riE 



£\,kfx,k — TLeo- 


(7b) 

(7c) 


In this Section we derive the macroscopic theory of 
electron transport in graphene in the absence of dissipa¬ 
tion. The resulting hydrodynamic equations represent a 
generalization of the Euler equation of an ideal liquid to 
Dirac fermions in graphene. 


Here we introduced the short-hand notation 


1. Continuity equations in graphene 


L 




Y. 


The hydrodynamic equations for the densities and cur¬ 
rents can be obtained by averaging the QKE ([^ with 
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respect to the modes (5d|. This yields the continuity 
equations for the hydrodynamic densities, 


9tn + V-j = 0, (11a) 


dtni + V • = 0, (11b) 


dtUE+ V ■ Je = eE ■ j, (11c) 


as well as the equation for the energy current 

“t" G-TiEq^ CTli^UXB^Q, — jE,al^dis’ (1^) 

Using the local distribution function Q, we can express 
the energy current in terms of the hydrodynamic velocity: 


3e — 


iUEU 
2 +m2’ 


(13) 


The equation (12) includes the momentum flux or 
stress tensor 


2. Hydrodynamics of ideal electron liquid 


In the traditional hydrodynamicthe ideal fluid is de¬ 
scribed by the Euler equation. The Euler equation is 
nothing but the continuity equation for the momentum 
density, where the stress tensor is expressed in terms of 
the velocity field. The latter is typically done on the basis 
of Galilean invariance. 

Similar equation can be formulated for the electron liq¬ 
uid in graphene. The momentum density is equivalent to 
the energy current which satisfies the continuity equa¬ 
tion (12). Substituting Eqs. (13) and (15) into Eq. (12) 


yields the Euler equation 


dn^Ma nEil-u"^) SuEU^up 
^* 0 I „,2 o I .,2 - + ^/9 o I „,2 (18) 


This equation is complemented by the continuity equa¬ 
tions ( [IT| ) and the self-consistency condition ( [T7| ). This 
set of equations generalizes the hydrodynamics of an ideal 
liquid to Dirac fermions in graphene in the absence of dis¬ 
sipation. 


n 


E 

OL^ 



£\,kVaVpf\,k ■ 


(14) 


C. Dissipative corrections 


In the absence of magnetic field, we use the distribution 
function (Q and ([^ to express 11®^ in terms of u: 

~ (15) 


Here the last term JH-® describes the dissipative effects 
that are considered in the next Section. The first term 
is the generalization of the usual stress tensor of an ideal 
liquicP to the case of Dirac fermions in graphene. The 
unusual form of Eq. (15) reflects the absence of Galilean 


as well as Lorentz invariance in the system. 

The electric and imbalance currents can be similarly 
related to the hydrodynamic velocity 


j = nu + 6j, 


(16a) 


jj = niu + 6jj. (16b) 

Here again we have introduced the dissipative correc¬ 
tions Sj, Sjj. Neglecting these terms along with 
the equations presented in this section describe the flow 
of the ideal electronic liquid. Since we are describing 
charged particles, the electric field should include the self- 
consistent electric (Vlasov) field 


In this Section we extend in the hydrodynamic the¬ 
ory of Dirac fermions in graphene by taking into account 
dissipative effects. We use the explicit form of the non¬ 
equilibrium distribution function (™ to evaluate the dis¬ 
sipative corrections 6 j, Sjj, and Miap- Comparing our 
results with the canonical form of the viscous terms in 
the stress tensor, we find the expression for the viscos¬ 
ity coefficients in graphene. We calculate the dissipative 
corrections to leading order in the gradient expansion. 
The parameter controlling the expansion is similar to the 
Knudsen number Kn = /hydro/^Vj where /v is the char¬ 
acteristic length scale of hydrodynamic fluctuations. 


1. Dissipative corrections to the currents 


Macroscopic equations that describe the electric and 
imbalance current densities j and jj can be obtained by 
integrating the kinetic equation similarly to the deriva¬ 
tion of Eq. (12). However, as j and jj are not conserved, 
the resulting equations contain non-vanishing contribu¬ 
tions of the collision integral. These contributions can be 
written in the form 


(m,/:/) = -(m,cj/(i)). 


(19a) 


r[r) =—'Vr Jcfr'V{r — r') Sn{r'). (17) 


iXv,Cf) = -iXv,CSf<^^^), 


where we have used a short-hand notation 

Here Sn{r) = n{r) — no is the local charge fluctuation, 
no is the background charge density, and V{r) = jr is 
the 3D Coulomb potential. 


( 5 )/)= / 9\,kfx,k- 

Jx.k 


(19b) 


(19c) 
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Using the distribution function (5b) we can now con¬ 


struct the explicit relation between the dissipative cor¬ 
rections to currents and the coefficients 



= M 


where the matrix JH is given by 


M = — 
2T 



(20a) 


vanishes (physically, because of the Galilean invariance). 
Technically, all macroscopic currents become equivalent 
and in particular are characterized by the same transport 
relaxation rate ^ which is much smaller than the 

usual rate ^ determining both the quasipar¬ 

ticle lifetime and thermalization. Further details of the 
calculation are relegated to Appendix |A 1| 

Solving Eq. (21) for the electric currents to leading 


order in the gradient expansion (i.e., in the Knudsen 


/ Cl 

Cx 

Ce/T\ 

Eqs. (16a) and (16b) 

c. 

Cl 

C\e\/T , 

(20b) 

VCe/T 

C\e\/T 

Ce'2 /T^ / 

( 


5j 

^Ji 




(24) 


+00 


Cx = NT J deiy{e)X • (20c) 


The coefficients C^/t, C|e|/t ^^nd 0^2/rp2 are proportional 
to the macroscopic densities 


a = 2n, 


q,| = 2n/, C ,2 = Zue/T. (20d) 


In Eq. (20c), T is the equilibrium background tempera¬ 
ture. 


The relation (20a) allows us to write the macroscopic 


equations for the electric and imbalance currents in the 
matrix form 


dt 


1 ( Vn — eEd^n 

2 \Vn/ - eEdfj,ni 


= -Cj 


Sj 
^31 


( 21 ) 


The matrix Cj plays the role of the collision integral in 
the reduced three-mode space. Its inverse is given by 


C~^ — 


Tl T2 
T3 Ti 


( 22 ) 


The transport scattering times Tj are obtained from the 
matrix elements of the linearized collision inte¬ 

gral C, where </> and ip' are the modes defined in Eq. (5d). 


The off-diagonal times r 2,3 change their sign for n —>■ —n. 
In the non-degenerate regime fj, T the times Tj are de¬ 
termined by temperature and electron-electron interac¬ 
tion, Tj = fj(fj,/T)/(agT), where fj{fi/T) is a smooth, 
dimensionless function. Close to the Dirac point. 


T2 = T3 = 0, 


(23a) 


while 


and 


T-^ = ^,,^—{v^,Cv^)^2.22alT, (23b) 


2r2 In 2 


where the vector i/j is given by 


1^.7 = 


^Vn^-iVn- 

iVn/- 


^ _ n 

Sue 

2ennj 
Sue 


- ^dp.ni 


(25) 


Here we have neglected the frequency dependence for¬ 
mally present in Eq. (21) since the hydrodynamic de¬ 


scription is valid at time scales much longer than the 
relaxation times due to electron-electron interaction that 


form the matrix (22). 


Individual terms in Eq. (25) allow for a simple physi¬ 


cal interpretation. The first term in each row describes 
the thermoelectric effect; the second term describes dif¬ 
fusion of electrons and quasiparticles; the last term leads 
to the finite conductivity of graphene due to electron 
interaction^^. The latter comprises a Drude-like term, 
which becomes more apparent if we identify the mass 
density p ^ 3nEl2n [see Eq. (41) and the text below] 
and a second term that gives rise to the finite conductiv¬ 
ity at the Dirac point for vanishing charge density n. 


2. Dissipative corrections to the energy stress tensor 


The macroscopic currents (10) are defined as the first- 


order moments of the distribution function with respect 
to the three modes (5d). The second-order moments yield 


the “generalized stress tensors” 


= / (l^lVaVpfxk- (26) 

Jxk 

Here the term with I = 3 is (up to the factor of T) the 
usual stress tensor ( [T4| ) . We also define the corresponding 
dissipative corrections 




r 1 

\sn, 

i = 1 

T 4 = 2T2 1n2^^^“’^^^“^ ^ 0.05 a^T. 


= / plVaVpSfxk = \ 


l = 2 (27) 

(23c) 

Jxk 


1 = 3 


Ear away from the Dirac point, /r ^ T, the system be¬ 
haves similarly to the usual Eermi liquid, where the trans¬ 
port mean-free time due to electron-electron interaction 


where the latter has been already defined in Eq. (15). 


The dissipative corrections (27) can be found by inte¬ 


grating the kinetic equation similarly to what was done 
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for the currents above. This way we find the relation 

(28) 


/TSUafi 

I Tsnip 





between and the coefficients from Eq. (5c). The 


matrix A4 is defined in Eq. (20b). Now we can express 


the right-hand side of the integrated kinetic equation in 
terms of the <511^*^. The resulting matrix equation reads 
[cf. Eqs. (|l^ and (21)] 


{(l^lVaVp.Cf) = -{(j}iVaVi3,Cdf^‘^'>) = -C^ . (29) 

Inverting the matrix collision integral we solve the 
above equation and find the dissipative corrections (27) 
similarly to Eq. (24): 




= cz 


^7r,a/31 


(30a) 


where to leading order in the gradient expansion 


^a /3 V • (nu) — Vanup — V i3nUa 
5a0V-{niu) - VaTliUp - V/jUjUa 
-{nEu) - VaflEUp - 


(30b) 

The matrix collision integral Ctt is discussed in detail in 
Appendix |A 2[ Hereafter, we restrict our discussion to 
the non-degenerate regime, /x <C T. Close to the Dirac 
point we find 


'C,r,ll 0 0 

= 2 I 0 C,r,22 C,r ,23 

0 C7r,32 C7r,33y 

with the matrix elements given by 

1 


C. 


■^5*^ — (AI ')kj- 


The traceless tensor is defined as 


(31a) 


(31b) 


(31c) 


Close to charge neutrality ( see A ppendix A 2 for details), 
all matrix elements in Eq. (31b) are of the same order 


1 




(32) 


— (MaPyCelap) ~ 


The dissipative correction to the stress tensor (15) is 
given by the third component of Eq. ( 30a| ) . To lead¬ 
ing order in the fluctuations of the densities, i.e. for 
Sue/ue ^ 1 as well as TStij/ue <SZ 1 and TSn/nE “C 1, 
the correction 6 Ii^ takes the canonical fornJS 

=-g\S/aUppUa-SanV-u], (33) 


with the viscosity coefficient 


r?=^(00l)C, 




" ) 

n/ 

V inE/2T I 


(34) 


see Eqs. (30). Close to the Dirac point this yields 


77 = T(T,r.in -f T^,2ni)/A + St^^sUe/S. (35) 


At the Dirac point the first term in Eq. ( |35[ ) drops out 
and we are left with two contributions to the viscosity g. 
The times r,r.i) T- 7 r .2 and T,r ,3 are obtained from invert¬ 


ing the collision integral (31) where the charge density is 
decoupled from the imbaL 


ance and energy densities: 


= 0 , 


E — ^ 1^71-,32 ^ 1 

2 C,r,23C7r,32 — C7r,22C7r,33 ’ 


1 


Ctt' 


2 C,r,22C 77,33 — C-n,2'3CTr ,32 
As a consequenc^AH 

r^{n = Q) = B T'^/c 
where the numerical coefficient is 


oc 


'5’ 


a^T 

ULg± 


5 oigTr-^a + 


9C(3) 

47r 




(36a) 


(36b) 


(36c) 


(37a) 


(37b) 


Here we have used the relations he = &C,{2>)T^/tt and 
nj = T^7r/3. Far away form t he D irac point we recover 
the usual Fermi-liquid viscositjPUl rj (xl /T^. 

Similarly to the classical hydrodynamicJi^, the viscos¬ 
ity is determined by the homogeneous equilibrium back¬ 
ground charge, imbalance and energy density, or equiva¬ 
lently by the equilibrium chemical potentials (/xo,±) and 
temperature T. The expression (33) implies vanishing 


bulk viscosity in graphene. This result is valid within 
the leading approximation in the virial expansion that 
justifies the kinetic equation ([^ as well as the distribu¬ 
tion function ([^. 


D. The canonical form of the hydrodynamic 
equations in graphene 


In this Section we combine the dissipative terms (33) 


and (24) with the equations of the ideal flow in graphene, 
see Sec. |IB 2[ The resulting theory generalizes the 
Navier-Stokes hydrodynamics to the Dirac fermions in 
graphene. 

The complete hydrodynamic description includes the 
equations of motion, continuity equations, and equa¬ 
tions of stat^. Within the local equilibrium approach 




























in graphene, the expression for the hydrodynamic pres¬ 
sure in terms of the energy density and the hydrodynamic 
velocity u is highly nonlinear 


(1 - u‘^)nE 

2 + u^ 


(38a) 


For small velocities the pressure assumes the standard 
value for a scale invariant gas, Pq = ueI'^, however, for 
large velocities approaching unity u ;< 1 it vanishes as 
~ (1 — u^). The enthalpy of the system VF = -f P is 
then given by 


momentum-dependent optical conductivity in graphene 
up to the subleading order in q/ur [and for l/(wTdis —t 0 )] 


cr(w, q) =ao + 


2 ie^n^ 

Sueuj 


Kp 


2 irjiu\ 

Sue 


)\ 


(41) 


iq^ r^2 
+ — 

UJ 


t'i + T 2 T 3 f 2 e^n^ 


Zue 




T 2 (ti + T 4 ) f 2 e^nni 

Sue 


+e^dnU 




w = 


2 w 

2-f u2’ 


w = ue + Pq = 3 n£;/ 2 , 


Here (Tq is the electron- electr on contribution to the dc 
(38b) conductivity in graphen d^^ l ^^ l 


with the latter being the linear enthalpy of graphene. 

The continuity equations 0 are now modified by the 
dissipative terms (25), 


dtn + V • (nit) = —V • 6 j, (39a) 

dtui + V • (n/it) = -V • Sjj, (39b) 

dtUE + V • (Wu) = enE ■ u. (39c) 


o-Q 



2n^ \ 


( d^ni 

V 2 



. (42a) 


In the above results, n, n/ and ue are the equilibrium 
background densities; the scattering times t, follow from 
Eq. ( [^ (see also Appendix |A2[ ). At the Dirac point, the 
electronic com pressibility in graphene is = 4Tln 2 / 7 r, 
and hencd^ 


(To = AeVckg, (42b) 


Finally, adding the dissipative part of the stress tensor 


(33) to the Euler equation (18) we obtain a generaliza¬ 


tion of the Navier-Stokes equation to Dirac fermions in 
graphene. Using the equations of state (38), we can bring 


the resulting equation to the canonical form (cf., Ref fTTll 
Wdtu + W{u ■ V)u + VP + udtP + u{ 6 j ■ E) (40) 

= en[E — u{u ■ E)] -|- rjS/^u. 


The term udtP in the left-hand side of Eq. (40) is re¬ 


flection of nearly relativistic nature of charge carriers in 
graphene. In the limit u —> 1, the electric field on the 
right-hand side of Eq. (40) does not affect the absolute 


value of the velocity which is limited by Vg. 

The complete system of the hydrodynamic equations 


in graphene includes Eqs. (38), (39), (40), as well as the 


equations defining the non-equili 
electric and imbalance currents (|24|). 


Drium corrections to the 


II. LINEAR RESPONSE 


A. Nonlocal optical conductivity 


where we find A = 0.19 (previously, the value A = 0.12 
was reported in Ref. | 


At (7 = 0, the conductivity (41) can be interpreted in 
terms of the usual Drude formula, where the role of the 
effective mass density is played by the ratio 3nE/{2n). 


The result (41) suggests a possibility to measure the 


viscosity coe fficien t in graphene in nonlocal transport 
measurement^^^^. However, precisely at the Dirac point 
(n = 0 ) the optical conductivity is independent of viscos¬ 
ity. Physically, viscosity is associated with the momen¬ 
tum density, i.e. the energy current. At the Dirac point, 
the energy and electric currents decoupl^^^ and hence the 
conductivity is unaffected by viscous effects. 

Finally, let us remark on the apparent contradiction 
between Eq. (41) and the corresponding result of Ref. |H1 


where it was found that the expansion of the optical con¬ 
ductivity in q/u! contains linear terms missing in Eq. (41). 
The reason for this disagreement is that we have cal¬ 
culated the response to the total electromagnetic field, 
while the result of Ref. |HI represents the response to the 
external field. In the latter case one has to take into ac¬ 
count screening which leads to the linear in q terms in 
nonlocal conductivity. 


Evaluation of the linear response transport coefficients 
within the hydrodynamic theory is straightforward. Lin¬ 
earizing the Navier-Stokes equation, we recover the linear 
response theory derived in Ref. [18] with the important ad¬ 
dition of time- and momentum-dependent contributions. 
Solving these equations, we find the expression for the 


B. Hydrodynamic energy waves and plasmons 


In a formally infinite system, the hydrodynamic theory 


(38) - (40) admits solutions in the form of collective en¬ 
ergy waves with the dispersion relation (which we obtain 











































9 



discussed for momenta that are large compared to the 
characteristic scales of both disorder and interaction. 

In a regular 2D electron systems, electric current is 
relaxed by disorder and as a result, the plasmon waves are 
damped at the lowest momenta. The plasmon dispersion 
is given bjB^ 

( * \ 1 2 

w wH-= -xqvp, 

\ '^dis / ^ 

such that for momenta smaller than the inverse Thomas- 
Fermi screening radius 


FIG. 1. (Color online) Energy wave dispersion ( |43| l, for dif¬ 
ferent chemical potentials. In the main panel (a), we compare 
the energy waves in an ideal fluid (dashed lines) to the dis¬ 
sipative (viscous) flow (solid lines). The inset (b) shows the 
effect of disorder scattering at small momenta. The curves 
are calculated for l/Truia = 0.001. 


as an expansion in q/w < 1 ) 

i , . agao . 2 f V , n + D 
-I- - iq -1-^- 


w(q) = - 


^Tiis 


\nE 




/ 4agn^ 

2 \ Snsq 


(43) 


These solutions can be interpreted as the hydrodynamic 
zero modes corresponding to poles in the response func¬ 
tions, see Appendix]^ Here we have also taken into ac¬ 
count weak disorder, which is absent in Eq. (40). 


For pure systems in the absence of dissipation the dis¬ 
persion relation ( [4^ greatly simplifies. At charge neu¬ 
trality (n = 0 ), the leading term is linear in q, 

w(n = 0,Tdis-t cxo,7?0) « Ugq/v^. (44) 

This acoustic energy wav^^^ is analogous to the long- 
wavelength oscillations in interacting systems of relativis¬ 
tic particleJS^ sometimes called “cosmic sound”. Suc h 
oscillations play an important role in astrophysic^^^^^. 

Away from charge neutrality, the collective modes of a 
pure system exhibit the square root spectrum typical for 
2D plasmons: 


w(Tdis oo,r/^0) 


2.agq 

Sue 


(45) 


Let us stress, that this mode is not the usual RPA 
plasmon. The crucial point is that the hydrodynamic 
description developed in this paper is valid at length 
scales much longer than the scale Zhydro, associated with 
electron-electron interaction, i.e. for very small mome nta 
q l^ydro- contrast, the usual RPA plasmonj ^^ l ^^ l are 




(46) 


As a result, for momenta much smaller than the inverse 
mean-free path the plasmon dispersion is purely imagi¬ 
nary, as expected for diffusive systems. For energy waves 
in graphene disorder scattering plays a similar role, see 
Eq. (@. 

Moreover, in graphene the electri c current is rela xed 

also by electron-electron interaction j^b8 | 24 | 36 | 37 | 48 | 4 9]^ 

a result, the plasmon modes are dampecP^ similarly to 
Eq. (46) even in the absence of disorder: 


UJ = — 


2 tp,, 


( 

O’O c^g D + 74 

1 ' VI 

2 



2e2 q 4 

2 rdisq^ / 


where 



where = >cq/2 for q x with the inverse Thomas- 
Fermi screening radius being x = 2TTaa(dan). Such plas¬ 
mons exist even at charge neutrahtjPlf (fQj- T > 0). Thus 
for small momenta, the plasmons are overdamped in con¬ 
trast to the energy waves (43). However, away from 


charge neutrality, the energy waves hybridize with the 
charge sector due to Vlasov self-consistency leading to 
dynamic oscillations of the charge density with the dis¬ 


persion (45), that is similar to ujp, but with a smaller 


prefactor. These oscillations should be expe rimen tally 
observable in the same way as usual plasmon^^Ulll^ pj-Q- 
vided that the samples (as well as the time scale of the 
measurements) are in the hydrodynamic regime. 

Far away from the Dirac point (/i 3 > T), the distinc¬ 
tion between the charge and energy sectors of the the¬ 
ory disappears, such that the energy waves coinci de w ith 
the usual plasmoiP^ for /r ^ T, the dispersion (45) re¬ 
produces ojp. Technically, the transport relaxation time 
due to electron-electron interaction that determines the 
above plasmon damping becomes much longer than the 
usual electron-electron scattering time that is responsible 
for thermalization in the system, see discussion following 
Eqs. (p3|. 


Viscous forces influence the collective modes (43) in 


the higher order in q/uj, cf. Eq. (41). Unlike the case of 
the optical conductivity, here viscosity enters in a linear 
combination with the scattering times ti and T 4 . Conse¬ 
quently, measuring the energy wave dispersion might not 
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be the best way to find the viscosity in graphene. How¬ 
ever, combining such measurements with the measure¬ 
ment of nonlocal conductivity, one can find experimental 
values for not only 77 , but also the scattering times 
The above results are illustrated in Figg where we 
plot the dispersion (43) for different chemical potential. 


The inset illustrates the role of disorder, cf. Eq. (46). 


III. NONLINEAR EFFECTS: RELAXATION OF 
A HOT SPOT 


ballistic propagation in an ideal system are determined 
by the initial conditions. 

The decay of the hot spot into the energy waves does 
not lead to an immediate relaxation of the initial energy 
density profile, see Fig. gb). In contrast to the three- 
dimensional flow, the Green’s function of the 2D wave 
equation exhibits a long-time tail, As a conse¬ 

quence the relaxation of the hot spot in the dissipation¬ 
less limit without Vlasov field shows power law decay. 
This slow relaxation of the energy density around the 
origin (afterglow) can be seen in Fig. ^b). 


In this Section we report results of a numerical inte¬ 
gration of the nonlinear hydrodynamic equations (38) - 
(40) describing relaxation of a hot spot. 


Let us prepare the system in a homogeneous, equilib¬ 
rium state characterized by the charge density (i.e.. 


(0) 


and 


away from charge neutrality), energy density 
imbalance density n. On top of this equilibrium back¬ 
ground, we create a hot spot: a locally elevated energy 
density. For simplicity, we choose a Gaussian profile with 
the peak height ue = 1.8n^\ see Fig. [^a). The result¬ 
ing non-equilibrium state will serve an initial condition 
for the subsequent time evolution that follows Eqs. (38) 
- (|40|. 


The computer simulations are performed in a semi- 
implicit schem^^. The diffusive and viscous corrections 
are discretized implicitly. This scheme is suitable for a 
wider class of problems that are characterized by compet¬ 
ing convective and diffusive terms. Moreover, the simu¬ 
lations are performed on a staggered grid to avoid un¬ 
physical density oscillation^^. 


A. Ideal flow 

We begin with the evolution of the hot spot in an ideal 
system described by the Euler hydrodynamics ([TT])-(1T§. 
Here we assume that the system is not subjected to any 
external fields. 


1. Pure energy flow 


2. Charge fluctuations 

In a charged system, i.e., in the presence of the self- 
consistent electric field, the cosmic sound wave shown 
in Fig. g is accompanied by fluctuations of the charge 
density, see Fig. g 

The excess energy density generates the pressure force 
described by V^Hj'^. This creates the initial energy flow 
that corresponds to the nonzero hydrodynamic velocity 
u, se e Eq. and hence translates into an electric cur¬ 
rent (I 6 a), which is coupled to t he ch arge density by 
means of the continuity equation ( |lla[ ). This way, the 
initial evolution of the excess energy density leads to a 
depletion of the charge density at the origin. 

Now, the non-equilibrium charge density profile results 
in the self-consistent electric field [due to Vlasov terms 
Remarkably, in the absence of dissipation the elec¬ 
tric field partially compensates the pressure force lead¬ 
ing to the appearance of a stable soliton-like composite 
density profile at the origin: after the initial outflow of 
energy carried away by the cosmic sound waves, some 
excess energy density remains at the point of the initial 
perturbation accompanied by the dynamically generated 
dip in the charge density, see Fig. g 

The establishing of the depletion in the charge density 
is accompanied by the charge flow shown in Fig. g Al¬ 
though that figure shows the flow in the presence of dis¬ 
sipation, at the short time scales used in the figure the 
dissipative effects are still weak and the resulting flow 
can be considered dissipationless. 


Within the hydrodynamic approach, the energy flow is 
coupled to the charge flow by means of the self-consistent 
electric field ( [T7| ). Turning off the Vlasov terms (i.e., 
setting E — 0), we arrive at an essentially neutral system 
where the energy flow is decoupled from the rest of the 
degrees of freedom. 

In such a system, creating an excess energy density 
leads to excitation of ballistic (due to absence of dissipa¬ 
tion) energy waves with the linear dispersion (44). This 
flow is illustrated in Fig.gb), where we plot the resulting 
energy density profile along the line y = 0 as a function 
of the ^-coordinate and time. In Fig. g we use arbitrary 
units, since the time and length scales associated with the 


B. Dissipative relaxation dynamics 

Gonsider the hot spot relaxation in a fully interacting 
system, i.e. in the presence of dissipation. We start with 
the same initial condition as before, but now the system 
evolves under the Navier-Stokes hydrodynamics (38) - 

pl). 


The hot spot evolution now proceeds in two stages. 
The first stage is similar to the ideal flow, where the 
quasi-stable charge-energy density profile is established 
at the origin. During this stage, some energy and charge 
are being carried away from the hot spot by the emit¬ 
ted energy waves. The metastable patterns, such as 
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FIG. 2. (Color online) Hot spot relaxation of a neutral ideal fluid obtained from the Euler hydrodynamics (111 - (181 without 
the Vlasov self-consistent electric field, E — 0. The left panel shows the initial energy bump with the height ue = 1.8n^\ The 
right panel shows the evolution of the energy density (in units of the equilibrium background, UE/n 
^-coordinate and time (arbitrary units) along the line y = 0. 


(o)^, 

E 


as a function of the 



FIG. 3. (Color online) A snapshot of the charge den¬ 
sity n. The equilibrium value of the charge density is 
= 1.9 X 10®cm“^. The initial height of the energy bump 
is n_B = 1.8n®. The inset (b) illustrates the soliton-like com¬ 
posite profile that is established at the origin. The blue curve 
shows the dip in the charge density and the red curve shows 
the excess energy density. The arrows show the balanced 
hydrodynamic forces: the pressure (red arrow) and the self- 
consistent electric field (blue arrow). 


the charge-energy complex in Fig. [^b) and the traveling 
waves are formed due to the nonlinear interplay between 
the charge and energy sectors. These patterns were sta¬ 
ble in the absence of dissipation, but now acquire a finite 
lifetime. 

Dissipative effects are characterized by a distinctly 
longer time scale compared to the initial evolution of the 
hot spot. These effects are manifested during the sec¬ 
ond stage of the hot spot evolution. Here the electron- 


Sn 



y=0 


0 

-8 -4 0 4 8 

X [/illl] 


FIG. 4. (Color online) The charge density as a function of x 
along the line y = Q for short enough time scales such that 
the system is effectively in the dissipationless limit. 


electron interaction leads to damping of the emitted 
waves, with the damping rate given by the imaginary part 
of the spectrum (43). In the clean limit, the dominant 
contribution to the damping rate is linear in q (similar to 
the 2D Maxwell relaxation, but with uo determined by 
electron-electron interaction). Furthermore, the soliton- 
like charge-energy complex is no longer stable and decays. 
However, the depletion of the charge density at the ori¬ 
gin remains visible for at least several picoseconds after 
the initial perturbation, see Fig. and hence shou ld be 
detectable by modern experimental technique^^HM] 
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IV. CONCLUSIONS 


In this paper we have presented a hydrodynamic de¬ 
scription of the electronic transport in graphene. Our for¬ 
malism allows for a consistent treatment of nonlinear hy¬ 
drodynamic effects as well as dissipative phenomena due 
to electron-electron interaction. Our theory describes the 
following hydrodynamic modes: the energy, particle and 
imbalance densities and the energy current. The electric 
and imbalance currents are relaxed by electron-electron 
scattering and have to be obtained from the equations of 
state. The resulting macroscopic description includes a 
generalization of the Navier-Stokes equation in graphene 
(40), the nonlinear relations (38) between the hydrody¬ 


namic pressure and enthalpy and the hydrodynamic ve¬ 
locity u that is related to the energy current. These 
relations play the role of the equations of state. Finally, 
the three macroscopic densities obey the set of continuity 
equations (39). 


Having derived the hydrodynamic theory from the 
Boltzmann kinetic equation, we are able to calculate ex¬ 
plicitly the set of scattering times that determine the 
coefficients in the hydrodynamic equations, in particu¬ 
lar the viscosity (37) and the dc-conductivity at charge 


neutrality (42b). The latter is the manifestation of the 


non-Galilean-invariant nature of the electronic system in 
graphene, where the electric current can be relaxed by 
electron-electron interaction. 

In laboratory experiments, viscous effects can be de¬ 
tected, for instan ce, b y measuring nonlocal conduc¬ 
tivity in graphene^SEil, Within linear response, vis¬ 
cosity affects the conductivity away from charge neu¬ 
trality and at nonzero momenta. Another experimen¬ 
tally detectable viscous effect is the plasmon lifetime in 
graphene. Although the viscosity coefficient enters the 
plasmon damping in a linear combination with other 
interaction-dependent parameters, see Eq. (43), measur¬ 


ing both the plasmon lifetime and nonlocal conductivity 
may give experimental access to several relaxation times 
determined by electron-electron interaction. 

Beyond linear response, we have considered the sim¬ 
plest example of nonlinear phenomena in graphene - the 
relaxation dynamics of a hot spot, see Fig. This anal¬ 
ysis takes into account the convective nonlinearities and 
the residual Coulomb interaction. In the macroscopic 
equations, the latter manifests the self-consistent electric 
field due to charge fluctuations and the dissipative cor¬ 
rections. We have found that the hot spot relaxation 
proceeds in two stages. The first stage, lasting no longer 
than few picoseconds, is characterized by the metastable 
charge-energy profile at the origin and the traveling en¬ 
ergy waves that carry excess energy and charge away 
from the hot spot, see Figs.|^|^ The emitted waves ex¬ 
hibit characteristic modulation due to the self-consistent 
Vlasov electric field. During the second stage, dissipative 
effects start playing a definitive role in the process leading 
to he diffusive charge propagation, damped energy waves, 
and the decay of the soliton-like charge-energy profile at 


the origin. The dissipative effects are much slower than 
the initial evolution of the hot spot. In particular, the 
metastable charge-energy profile remains visible at times 
of order lOps, which should be detectable in laboratory, 
see Fig. 

The traveling energy waves are accompanied by fluc¬ 
tuations of the charge density due to nonlinear coupling 
between the energy and charge sectors in the theory away 
from charge neutrality. Precisely at the Dirac point, the 
energy waves have linear dispersion (44), similar to the 
cosmic souncP^. For finite background charge densities 
the dispersion of the energy waves (45) becomes similar 
to the usual 2D plasmons^Sl 


with its intrinsic life-time 
determined by electron-electron interaction. However, as 
the hydrodynamic theory is valid only for time and length 
scales that are much larger than the typical scales associ¬ 
ated with the electron-electron scattering, the true plas¬ 
mon modes remain overdampecP^. However, far away 
from charge neutrality ^ T) we recover the usual plas¬ 
mon in graphene. 

The hydrodynamic theory presented in this paper is 
valid as long as quasiparticle recombination processes re¬ 
main slow (technically, infinitely slow). At time scales 
exceeding the recombination times the imbalance density 
is no longer conserved and the structure of the hydrody¬ 
namic equations changes. However, the Navier-Stokes 
equation (40) is independent of the imbalance density 
and remains valid even at the longest time scales. 

The problem of the hot spot relaxation and trav¬ 
eling energy waves considered in this paper is closely 
related t o recen t experimental imaging of plasmons in 
graphen^2ll2lI131 while the existing experiments are fo¬ 
cusing on the high-frequency optical phenomena, we hope 
that our investigation of the energy waves in graphene 
will motivate future measurements in the low-frequency, 
hydrodynamic regi me. A t the same time, nonlocal trans¬ 
port measurementd^oiMJ nray uncover exciting manifesta¬ 
tions of the nonlinear, viscous flow in graphene including 
vortices and laminar wake. 


Our hydrodynamic theory can be further applied to 
more realistic, experimentally relevant geometries in or¬ 
der to study possible realizations of the plethora of hy¬ 
drodynamic phenomena in graphene. After a straightfor¬ 
ward generalization, the theory allows us to consider the 
thermoelectric effects as well as the effects of the external 
magnetic field. This work will be reported elsewhere. 
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Appendix A: The ee-collision integral 

The electron-electron collision integral in the QKE (§ is given by 

Stee[f]= ^ f \M\‘^{2TT)^S{exp + ei,k-ex'p'-e^>k>)S{p + k-p'-k') (Ala) 

X {fx'.p'U'M [1 - f,y,p] [1 - fx.k] - fx,pfv,k [1 - fv',p'] [1 - f^',k'] } ■ 

Here the matrix element of Coulomb scattering is given by 

|M|2 = N\Viuj,q)\‘^exp,yp'e.k,.'k', (Alb) 

with the graphene specific Dirac factors 


0A,p;A',p' — 2 pp' ) — 2 ■ 


(Ale) 


prohibiting backscattering. In Eq.(Alb), uj = ex,p — ^x'y is the transfered energy and q = p' — p - the transfered 
momentum. 

Linearizing the collision integral ( |A1| wit h respect to the deviations of the distribution function from the local 
equilibrium Q, we obtain the operatoi*211SI 

CSfx.k= X! f \Mf{2TTfS{ep + ek-ep'-Sk')S{p + k-p'-k') (A2) 

Jk,p',k' 

[Sfx,p + SU^k - Sfx',p' - ■ 


^JX,pJv,k 


1 - 

JA'.p' 


1 - 

v' ,k 


1. Transport scattering times due to electron-electron interaction 


In this section we give explicit expressions for the scattering times Ti constituting the matrix collision integral in 
the space of macroscopic currents j and jj, see Eq. (21). These equations are obtained by averaging the QKE with 
respect to v and An. Therefore the right-hand side of Eq. is given by 


Cj 


Sj\ _ /(n,C5/W)A 
SjjJ {{Xv,C6f(^^)J 


(A3) 


The scalar product (•, •) was defined in Eq. (19c). 


Dissipative corrections to the macroscopic currents 6 j and Sj j are determined by the non-equilibrium contribution 
to the distribution function ([^ as 


5jk = {cjxkV,-T5f^^^d,f^^'^), 


(A4) 


such that 


Sj = Sj^, Sjp = Sj 


2- 


Here we remind the reader that the terms proportional to u in Eqs. (161 follow directly from the local equilibrium 
distribution 0- The functions are the modes ( |5d[ ). 

Now we can use the definition (A4| to express the coefficients non-equilibrium distribution (Im in terms 

of Sj and Sjj. This allows us to find the explicit form of the matrix collision integral Cj, see Eq. (22). Mter some 
algebra, we find 


2 

^ ^ [A4 ]j7c (^/Pa: Pa) 1 


(A5) 


where the matrix A4 is given by Eq. (20b I. 
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r(o)r(2) 




FIG. 5. a) The Aslamazov-Larkin-type diagram corresponding to the term The product comprises the 

Maki-Thompson-type diagram b) as well as self-energy correction c). 


The matrix elements in Eq. (A51 can be evaluated explicitly using the methods of Refs. [T8| and [24l Noting that in the 
integrated electron-electron collision integral the summation over scattering states {|A, fc), |A', fc')} and {\v,p), \v',p')} 
separates, we express the matrix elements as 


1 


{(j)Va,C(l)'vp) = — I duj I (fq 


sinh2(w/2T) L 


r00qa/3(w, (w, q) - q)Vyj^{u:, q) 


<f>,a V 


Here the vertex functions are defined as [A' = sign(e;^^p -|- w)], 

= ^ f S{ex,p - ey^p+q+u) (/® - fy]p+q) 0Ap;A'.p+q, 

^ JX.v ^ ' 


(A6) 


(A7a) 


— rp <5 (£a,p ^A',p+q + w) p /Icp-t-q) ®^.pA'.p+q [^^'.p+q^^Ap ^Ap] 


A2) 

00' 


•',ap~rp f ^A',p-|-q + w) p /A',p-|-q) ®A,p;A',p+q 

^ 'J \,P ^ ' 


(A7b) 


(A7c) 


X [<^A',p+q 'CACp-l-q 4>Xp ^’Ap] ^ \4>y ^p+q 'V\',p+q 4>Xp 'I'Ap] 

The product can be represented with the help of the Aslamazov-Larkin-type diagram in the Boltzmann limit, 

whereas the product contains the Maki-Thompson-type diagrams as well as self-energy corrections, see Fig.j^ 


The resulting values ( |23[ ) are most conveniently calculated in the local co-moving frame, where the hydrodynamic 
velocity entering the local equilibrium distribution functions in Eqs. (A7) vanishes. The obtained results are then 


valid in arbitrary reference frame based on the principle that the relaxation are independent of the reference frame 
(generalizing the Galilean invariance to the arbitrary spectrum). 


2. Dissipative corrections to the stress tensor 


The collision integral C,r can be calculated along the same lines as Cj in the previous Section. Averaging the QKE 
with respect to the tensor quantities such as we find the contribution of the collision integral in the form similar 

to Eq. (A31 


/ \ 
Ctt dlii^aP 
\T-HIIe,o.p) 


The stress tensors were defined in Eqs. (26) and (27). 
Defining the deviations from equilibrium as 


{VaVp.CSf^'^'^) \ 

{XVaVp,CSf^‘^'') 

isvo.vp/T,C6f(^^) J 




(AS) 


(A9) 
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we can express the coefficients in the non-equilibrium distribution function (j^ in terms of . Similarly to the 
arguments presented in the previous Section, this yields the explicit form of the matrix collision integral 


[C^7r]/fc — ^ ^ ^ [.Ad ]j/c H■ 
1 = 1 


(AlO) 


Here the matrix M. is given by Eq. (20b I and the traceless tensor is defined in Eq. (31c). 

The matrix elements {(f>iIap,C(f>jIap) can be evaluated similarly to Eqs. ( |A7[ ): 

Here the tensor vertex functions are [A' = sign(eA,p — w)], 

~ T Jx ~ ^"^bp-l-q + (/a,p ~ /Iqp+q) ®A,p;A',p+q [<py,p+q IaP,p+q ~ 4’\p ^a/3,p] ) 


7 ( 2 ) 

"“(pep' ,ct0'y5 


{oj,q) — ^ J S{e\^p sy ,p+q + oj) (^fx p /I'.p+q) ®->'.pA'>p+q 


(All) 

(Al2a) 

(A12b) 


X ['?^’A',p+q IaP,p+q 4'Xp IaP,p\ [<(’A',p4-q ^'lS,P+q 4'Xp -^7<5.p] ' 

For further calculations it is useful to express the tensor lap in terms of the basis vectors {q = q/q,q± = z x q}, 

^a/3 ^k.qi^QaQ^ ^ajs') “I” (Q_L,aQ/3 “1“ 9aQ_L,/3) 5 

where 

{k • q^){k • q) 


^k.q — 


(^■ 9 )" iV 1 1 I 1 


ddk.q — 




(A14) 


Due to the conservation laws of the electron-electron interaction we effectively have A —> A. Using the d-function in 
Eqs. (A12al and (A12b) one obtains (e = eA,fc), 


Ak.q = (o;" - q^) 


2 2N(2£ + w) 2 -g 2 


8e^g2 


(A15) 


Furthermore, the coefficient B drops out in the vertex function since it is antisymmetric in the angle between 


( 2 ) 

and k. In the tensor vertex function ^ap-ysi^^’ d) g^t ^ separate contribution from A and B but they are orthogonal. 
For B we obtain with the help of the d-functions {e = £>,*)) 

x/ {q^ — [(2s -I- ujY — q^] (w^ — q^ — 2ew) 


Bk.q = sign(fe • gj_)- 


4^2q,2 


(A16) 


Finally, with the help of the angular averages 

dipq qa^pq^QS — ^(^a7^/35 4“ ^aS^pj Sapdjs') •> 


j d(pq (^qaqp daP^i^q^qX — ^{dajdps ~h dasdpj dapdjs}: 


d(pq {qj_^aqp 4“ qaqj-,p^^qA.,xQs 4- q^Q-L^s') — ^{da^^ps 4- dasdp^ dapd^s^ 


and the projected quantities obtained after averaging Eqs. (A12) over the angle ipq of the transfered momentum 

q, 


^(1) 


{u!,q) — J' J ^-'bp+q 4 ”(/a,p fx'^p+q^ ^x,p-,x',p+q[<Px',p+q A^+q.q (f’XpAk^qP (Al7a) 
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"( 2 ) 


/(w,g) 


1 

T 


f - 

J A,p 

X \4’\'-p+q 


■ ^A'.p+q + w) (^fxp f\’]p+q^ 0A,p;A',p+q 

Ap+q,q — 4>\p Ap^q] ['5^A',p+q ^p+q,q ~ ^Xp Ap,q] ) 


(A17b) 


"( 2 ) 

-_L,00 





SA'.p+q + ‘^) (/a,p 


AO) 

J X' ,p+q 


Qa,p;A',p+<7 


(A17c) 


X \4'X' ,p+q Bpj^q q (pXp 7?p,q] [(/)^/ ,p+q Bp^q q (j)^p Bp.q^ , 

we can write the matrix elements as 

{(j)Ia0,C4>'ljs) = —^{Sa^Sps + SasSp'Y — SapSjs) [duj [(f q } ^ (A 18 ) 

Iott j j smh {ijj/21 ) 

(^> (w, g) + (w, g)r(°) (w, g) - (w, g)5^V (w, g) . 

Here we can drop the terms proportional to Sajs since the energy stress tensors are traceless. Due to their symmetry 
in a o /? we effectively have 




(A19) 


The matrix elements (A18) determine the quantities CT^^ij, Eqs. which in turn determine the viscosity (|35[). 


Appendix B: Linear response functions 


In linear response we linearize the hydrodynamic equations with respect to the linear fluctuations of the hydrody¬ 
namic quantities, 6n, 5 ni, Sue, Su\ 


n^n + 5n, ni^nj + 5ni, ue ^ nE + Sue, u ^ Su. 
We furthermore introduce the response functions to the external perturbation E = —iqqx 

6n = XnA, Sni = xiA, Sue = TxeA, Su = -i-^^XuA- 

quE 


(Bl) 


(B2) 


Linearizing the continuity equations (391 and the Navier-Stokes equation (401, we find the matrix equation for the 
response functions Xi'- 


(-ICU + ^g2 - 27regao ^g^ - (q^ 


V 


'yh — STrego-Q 

0 

_ 47re^n 

3T ^ 


-iuj -f ^g2 - 
0 
0 


3nE 

n/T4+nr3 2 

Sue • ^ 
—iuj 
_q 

3 


-iuj -f T^. ^ -b 


—g 

riE ^ 
niT 
riE ^ 

h 
-1 
Mis 







-g^CTo/e 

-q^olje 


—a 

3T y 


(B3) 


where [cf. Eq. (42a)] 


(To = e 


n 


dpU 


2n^ 

Sue 


+ T2 


/ dpUj 2nni \ 

V 2 Sue J 


T3 


dpU 


211? \ 
3nE, 


\ +74 


dpTii 


2nni \ 

Sue ) 


The dispersion (43) of the collective modes follows from zeros of the determinant of the matrix in the left-hand side 
ofEq. dB^ . 

In contrast to the energy waves and plasmons, which describe the response of the system to an external perturbation, 
the conductivity of an infinite system is defined as the response to the total electric field. Consequently, in order to 
find the conductivity (41), we need to consider the irreducible response functions, which satisfy the equation similar to 
Eq. (B3|, but without the Vlasov terms in the left column of the matrix in the left hand side. Then the conductivity 

nT 


is found from the Ohm’s law 

Sj = 


/ , 1 / , . nTi+ niT2 

0'0/e -b -(tiX„ -b T 2 X 1 ) -T- Txe +- Xn 

2 SriE qriE 


(B4) 
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where ip is now the total potential in the system (including the self-consistent Vlasov contribution). 


^ L. D. Landau and E. M. Lifshitz, Fluid mechanics 
(Butterworth-Heinemann, 2000). 

^ E. M. Lifshitz and L. P. Pitaevskii, Physical Kinetics 
(Butterworth-Heinemann, 1981). 

® D. Vollhardt and P. Wolfle, The superfluid phases of He¬ 
lium 3 (Taylor and Francis, 1990). 

^ K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, 
Y. Zhang, S. V. Dubonos, 1. V. Grigorieva, and A. A. 
Firsov, Science 306, 666 (2004). 

® S. A. Hartnoll, P. K. Kovtun, M. Muller, and S. Sachdev, 
Phys. Rev. B 76, 144502 (2007). 

® L. Fritz, J. Schmalian, M. Muller, and S. Sachdev, Phys. 
Rev. B 78, 085416 (2008). 

^ M. Miiller, L. Fritz, and S. Sachdev, Phys. Rev. B 78, 
115406 (2008). 

® M. Muller and S. Sachdev, Phys. Rev. B 78, 115419 (2008). 
® M. Muller, L. Fritz, S. Sachdev, and J. Schmalian, AIP 
Conference Proceedings 1134, 170 (2009). 

M. S. Foster and 1. L. Aleiner, Phys. Rev. B 79, 085415 
(2009). 

M. Muller, J. Schmalian, and L. Fritz, Phys. Rev. Lett. 
103, 025301 (2009). 

D. Svintsov, V. Vyurkov, S. Yurchenko, T. Otsuji, and 
V. Ryzhii, Journal of Applied Physics 111, 083715 (2012). 
A. Tomadin and M. Polini, Phys. Rev. B 88, 205426 
(2013). 

D. Svintsov, V. Vyurkov, V. Ryzhii, and T. Otsuji, Phys. 
Rev. B 88, 245444 (2013). 

T. V. Phan, J. C. W. Song, and L. S. Levitov (2013), 
arXiv: 1306.4972 (unpublished). 

L. S. Levitov, A. V. Shtyk, and M. V. Feigelman, Phys. 
Rev. B 88, 235403 (2013). 

A. Tomadin, G. Vignale, and M. Polini, Phys. Rev. Lett. 
113, 235901 (2014). 

B. N. Narozhny, 1. V. Gornyi, M. Titov, M. Schiitt, and 
A. D. Mirlin, Phys. Rev. B 91, 035414 (2015). 

A. Principi, G. Vignale, M. Carrega, and M. Polini (2015), 
arXiv: 1506.06030 (unpublished). 

S. A. Maier, Nature Phys. 8, 581 (2012). 

A. N. Grigorenko, M. Polini, and K. S. Novoselov, Nat 
Photon 6, 749 (2012). 

K. Damle and S. Sachdev, Phys. Rev. B 56, 8714 (1997). 
A. Principi, G. Vignale, M. Carrega, and M. Polini, Phys. 
Rev. B 88, 195405 (2013). 

M. Schiitt, P. M. Ostrovsky, M. Titov, 1. V. Gornyi, B. N. 
Narozhny, and A. D. Mirlin, Phys. Rev. Lett. 110, 026601 
(2013). 

J. C. W. Song and L. S. Levitov, Phys. Rev. Lett. Ill, 
126601 (2013). 

M. Titov, R. V. Gorbachev, B. N. Narozhny, T. Tu- 
dorovskiy, M. Schiitt, P. M. Ostrovsky, 1. V. Gornyi, A. D. 
Mirlin, M. 1. Katsnelson, K. S. Novoselov, et ah, Phys. 
Rev. Lett. Ill, 166601 (2013). 

D. N. Basov, M. M. Fogler, A. Lanzara, F. Wang, and 
Y. Zhang, Rev. Mod. Phys. 86, 959 (2014). 


Z. Fei, A. S. Rodin, G. O. Andreev, W. Bao, A. S. 
McLeod, M. Wagner, L. M. Zhang, Z. Zhao, M. Thiemens, 
G. Dominguez, et ah. Nature 487, 82 (2012). 

J. Chen, M. Badioli, P. Alonso-Gonzalez, S. Thongrat- 
tanasiri, F. Huth, J. Osmond, M. Spasenovic, A. Centeno, 
A. Pesquera, P. Godignon, et ah. Nature 487, 77 (2012). 
D. A. Abanin, S. V. Morozov, L. A. Ponomarenko, R. V. 
Gorbachev, A. S. Mayorov, M. 1. Katsnelson, K. Watan- 
abe, T. Taniguchi, K. S. Novoselov, L. S. Levitov, et ah. 
Science 332, 328 (2011). 

T. Taychatanapat, K. Watanabe, T. Taniguchi, and 

P. Jarillo-Herrero, Nature Phys. 9, 225 (2013). 

P. Alonso-Gonzlez, A. Y. Nikitin, F. Golmar, A. Centeno, 
A. Pesquera, S. Vlez, J. Chen, G. Navickaite, F. Koppens, 
A. Zurutuza, et ah, Science 344, 1369 (2014). 

A. Woessner, M. B. Lundeberg, Y. Gao, A. Prin¬ 
cipi, P. Alonso-Gonzalez, M. Carrega, K. Watanabe, 
T. Taniguchi, G. Vignale, M. Polini, et ah. Nature Ma¬ 
terials 14, 421 (2015). 

L. A. Ponomarenko, A. K. Geim, A. A. Zhukov, R. Jalil, 
S. V. Morozov, K. S. Novoselov, 1. V. Grigorieva, E. H. 
Hill, V. V. Cheianov, V. 1. Falko, et ah. Nature Phys. 7, 
958 (2011). 

R. Decker, Y. Wang, V. W. Brar, W. Regan, H.-Z. Tsai, 

Q. Wu, W. Gannett, A. Zettl, and M. F. Crommie, Nano 
Letters 11, 2291 (2011). 

A. B. Kashuba, Phys. Rev. B 78, 085415 (2008). 

M. Schiitt, P. M. Ostrovsky, 1. V. Gornyi, and A. D. Mirlin, 
Phys. Rev. B 83, 155441 (2011). 

A. Tomadin, D. Brida, G. Cerullo, A. C. Ferrari, and 
M. Polini, Phys. Rev. B 88, 035430 (2013). 

A. A. Kozikov, A. K. Savchenko, B. N. Narozhny, and A. V. 
Shytov, Phys. Rev. B 82, 075424 (2010). 

D. E. Sheehy and J. Schmalian, Phys. Rev. Lett. 99, 
226803 (2007). 

S. Chapman, Phil. Trans. R. Soc. Lond. A 216, 279 (1916). 

S. Chapman, Phil. Trans. R. Soc. Lond. A 217, 115 (1918). 

D. Enskog, Arkiv Mat. Astr. Fys. 16, 60 (1921). 

P. Arnold, G. D. Moore, and L. G. Yaffe, Journal of High 
Energy Physics 2000, 001 (2000). 

R. Sunyaev and Y. B. Zeldovich, Astrophysics and Space 
Science 7 (1970). 

P. J. E. Peebles and J. T. Yu, The Astrophysical Journal 
162 (1970). 

G. Zala, B. N. Narozhny, and 1. L. Aleiner, Phys. Rev. B 
64, 214204 (2001). 

E. H. Hwang and S. Das Sarma, Phys. Rev. B 75, 205418 
(2007). 

T. Stauber, J. Phys.: Condens. Matter 26, 123201 (2014). 

U. M. Ascher, S. J. Ruuth, and B. T. R. Wetton, SIAM 
Journal on Numerical Analysis 32, 797 (1995). 

M. Griebel, T. Dornsheifer, and T. Neunhoeffer, Numerical 
Simulation in Fluid Dynamics (Society for Industrial and 
Applied Mathematics, 1997). 



